% This file take the steady-state the file "todynaremone.mat"
% It writes the dynare file "code_dynare_rules.mod" which is the dynamic model
% to use perturbation techniques, using fiscal rules
% it doubles check that the steady state is OK, thanks to the "resid" function

% 


clear
warning('off','all');
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

isOctave && warning('off', 'Octave:array-as-logical');

load ../toDynare

trunc  = Ramsey.truncatedModel.truncatedAllocation;
lagMul = Ramsey.lagrangeMult;

Ntot  = Ramsey.truncatedModel.Ntot;    %total number of bins
Pi_h  = trunc.Pi_h';
S_h   = trunc.S_h;
K     = solution.K(1);
Ltot  = solution.L(1);
alpha = eco.alpha(1);
delta = eco.delta(1);
Y     = K^alpha*Ltot^(1-alpha);
FL    = (1-alpha)*K^alpha*Ltot^(-alpha);
xsis  = Ramsey.truncatedModel.xsis;

N1 = max(find( cumsum(S_h)<=1/3)); % type 1 history
N2 = max(find( cumsum(S_h)<=2/3)); % type 2 history

%xsis.xsiuE=xsis.xsiu1*0+1;
%xsis.xsiuE=xsis.xsiu1;


C1 = sum(S_h(1:N1).*trunc.c_h(1:N1));
C2 = sum(S_h(N1+1:N2).*trunc.c_h(N1+1:N2));
C3 = sum(S_h(N2+1:Ntot).*trunc.c_h(N2+1:Ntot));

coefB = 0.1;


ys    = eco.ys;


rho_u = 0.1;
G0 = .01*(1-rho_u); 


l     = Ramsey.l;
taut =  (( 1.0/eco.phi + 1.0)*eco.tau)/( 1.0/eco.phi +1.0 - eco.tau);
x = trunc.c_h - (1/(eco.chi*(1/eco.phi + 1)))*(trunc.l_h).^(1+1/eco.phi);


P     = ones(Ntot,1);
P(trunc.ind_cc_h) = 0;


file = 'code_dynare_trunc_rules.mod';
fid  = fopen(file, 'w');
formatSpec = '%25.20f';



%Weight = planner.weightp;%weigth_est;
omega_h = Ramsey.weights.omega_h;     %per history
omegab_h = Ramsey.weights.omegab_h;   %per history
%~ omega_h = Ramsey.weights.omega_param_h; #per history
%~ omegab_h = Ramsey.weights.omegab_param_h; #per history
%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     variables and params     %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['var\n'] ; fprintf(fid, str);

for h = 1:Ntot
    str = ['a',num2str(h),' '] ; fprintf(fid,str);     %  end-of-period
    if mod(h,10)==0; %multiple de 10
        str = ['\n'] ;         fprintf(fid, str);
    end;

end;


str = ['u r w FK FL A B K TFP Ltot Ctot C1 C2 C3 C1t C2t C3t Y I It\n'] ; fprintf(fid, str);
str = ['ut Kt Ct Bt wt rt taut tauk Zt Yt Lt tkt kappa BoYt G l T\n'] ; fprintf(fid, str);

%str = ['r w   mu   \n'] ; fprintf(fid, str);
%str = ['tau     \n'] ; fprintf(fid, str);


for h = 1:Ntot
    %str = ['a',num2str(h),' '] ; fprintf(fid,str);     %  end-of-period
    str = ['at',num2str(h),' '] ; fprintf(fid, str);   % beginning-of-period
   % str = ['l',num2str(h),' '] ; fprintf(fid, str);   % beginning-of-period
    str = ['x',num2str(h),' '] ; fprintf(fid, str);
    %str = ['psib',num2str(h),' '] ; fprintf(fid, str);
    %str = ['lambdacb',num2str(h),' '] ; fprintf(fid, str);
    %~ str = ['lambdac'e,num2str(h),' '] ; fprintf(fid, str);
    %~ str = ['lambdal',num2str(p),' '] ; fprintf(fid, str);
    %str = ['lambdacbt',num2str(h),' '] ; fprintf(fid, str);    % last period
    %str = ['lambdalb',num2str(h),' '] ; fprintf(fid, str);    % last period
    if mod(h,10)==0; %multiple de 10
        str = ['\n'] ;         fprintf(fid, str);
    end;

end;

str = ['; \n\n'] ; fprintf(fid, str);
str = ['varexo eps; \n\n'] ; fprintf(fid, str);
str = ['parameters\n'] ; fprintf(fid, str);
str = ['beta alpha phi abar delta gamma chi rho_u Gss;\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        initialisation        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['alpha','   = ',num2str(alpha,formatSpec),';\n']; fprintf(fid, str);
str = ['beta','    = ',num2str(eco.beta,formatSpec),';\n']; fprintf(fid, str);
str = ['phi','     = ',num2str(eco.phi,formatSpec),';\n']; fprintf(fid, str);
str = ['abar','    = ',num2str(eco.a_min,formatSpec),';\n']; fprintf(fid, str);
str = ['delta','   = ',num2str(eco.delta,formatSpec),';\n']; fprintf(fid, str);
str = ['gamma','   = ',num2str(eco.sigma,formatSpec),';\n']; fprintf(fid, str);
str = ['rho_u','   = ',num2str(rho_u,formatSpec),';\n']; fprintf(fid, str); %HERE
str = ['Gss','     = ',num2str(solution.G,formatSpec),';\n']; fprintf(fid, str);
str = ['chi','     = ',num2str(eco.chi,formatSpec),';\n']; fprintf(fid, str);

%%
equa = 1; % Numbering equations
tol = 0*10^-8;  % threshold for transition, to avoid to consider very small transitions

str = ['\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       model definition       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['model;\n\n']; fprintf(fid, str);
for h = 1:Ntot
    
    %%%%%%%% Budget constraint
    str = ['x',num2str(h),' = (1/(chi*taut))*l^(1+1/phi)*(',num2str(trunc.y0_h(h),formatSpec),')^taut  + T + (1+r)*','at',num2str(h),'-a', num2str(h), ...
           '; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    %%%%%%%% Definition of atilde (at)
    str = ['at',num2str(h),' = 0 ']; fprintf(fid, str);
    for hi = 1:Ntot
        if (Pi_h(h,hi))>tol   % hi to h
          str = ['+ ',num2str(S_h(hi)*Pi_h(h,hi)/S_h(h),formatSpec),'*a',num2str(hi),'(-1)']; fprintf(fid, str);
        end
    end
    str = ['; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    %%%%%%%% Euler equation
    if P(h)==1 % not constrained
        str = [num2str(xsis.xsiuE(h),formatSpec),'*(x',num2str(h),')^-gamma =', num2str(trunc.resid_E_h(h),formatSpec),' + beta*(1+r(+1))*(']; fprintf(fid, str);
        for hp = 1:Ntot
            if (Pi_h(hp,h)*xsis.xsiuE(hp))>tol
                str = ['+ ',num2str(Pi_h(hp,h)*xsis.xsiuE(hp),formatSpec),'*(x',num2str(hp),'(+1) )^-gamma']; fprintf(fid, str);
            end
        end
        str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    else
        str = ['a',num2str(h),' =  ',num2str(trunc.a_end_h(h),formatSpec),';']; fprintf(fid, str);str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    end

    
   
   
end

str = ['\n\n %% FOC planner \n\n']; fprintf(fid, str);



%%%%%%%%%%%%%%%% Aggregate Constraints %%%%%%%%%%%%%%%%

%%%%%%%% Resource contraint
str = ['G = Gss*(1+u);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['I = K - (1-delta)*K(-1);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['G  + Ctot + I = Y;\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%str = ['G + (1+r)*B(-1) + (Ctot +A - (1+r)*A(-1))  = Y +B - (r+delta)* K(-1);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Financial market clearing
str = ['A = B + K; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Production side
str = ['FL = (1 - alpha)*TFP*(K(-1)/Ltot)^alpha; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['FK = alpha*TFP*(K(-1)/Ltot)^(alpha-1)-delta; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['u = rho_u*u(-1) + eps;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['TFP = 1; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Y = TFP*K(-1)^alpha*Ltot^(1-alpha);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate labor supply
%str = ['Ltot = l*('] ; fprintf(fid, str);
%for h = 1:Ntot
%    str = ['+', num2str(S_h(h),formatSpec),'*',num2str(trunc.y0_h(h),formatSpec),'^(1+ taut/(1/phi+1))'] ; fprintf(fid, str);
%    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
%end;
%str = [');   %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;



%%%%%%%% Aggregate labor supply
str = ['Ltot = l*('] ; fprintf(fid, str);
for j = 1:eco.ny
    str = ['+ ', num2str(eco.Sy(j),formatSpec),'*',num2str(eco.ys(j),formatSpec),'^(1+ taut/(1/phi+1))'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [');   %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;



%%%%%%%% Aggregate consumption
str = ['Ctot = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+',num2str(S_h(h),formatSpec),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(trunc.y0_h(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate consumption
str = ['C1 = 0'] ; fprintf(fid, str);
for h = 1:N1
    str = ['+',num2str(S_h(h),formatSpec),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(trunc.y0_h(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

str = ['C2 = 0'] ; fprintf(fid, str);
for h = (N1+1):N2
    str = ['+',num2str(S_h(h),formatSpec),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(trunc.y0_h(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


str = ['C3 = 0'] ; fprintf(fid, str);
for h = (N2+1):Ntot
    str = ['+',num2str(S_h(h),formatSpec),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(trunc.y0_h(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;





%%%%%%%% Wage
str = ['w =  (1/chi)*(1/taut +1/(1/phi+1) ) *l^( (1+1/phi)*(1+1/phi)/(1/phi+1 +taut) )  ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Taxes
str = ['tauk = 1-r/FK;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['kappa = w/(FL^( (1/phi+1)*taut/(1/phi+1+taut) ));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate savings
str = ['A = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+', num2str(S_h(h),formatSpec),'*a',num2str(h)] ; fprintf(fid, str);
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%%% Fiscal system

str = ['tauk  = steady_state(tauk);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
%str = ['kappa = steady_state(kappa);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['taut = ',num2str(taut,formatSpec),';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

str = ['kappa = steady_state(kappa) -',num2str(coefB,formatSpec),'*(B - steady_state(B)) ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%str = ['T = -',num2str(coefB,formatSpec),'*(B - steady_state(B)) ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['T = 0 ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;



%%%%%%%%%%%%%%%% log-deviation variables %%%%%%%%%%%%%%%%
% There are denoted with 't': for variable x, xt = (x -x_ss)/x_ss, where x_ss is the steady state value of x.

str = ['ut = 100*u;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Kt = 100*(K/',num2str(solution.K(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Ct = 100*(Ctot/',num2str(solution.C(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['C1t = 100*(C1/',num2str(C1,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['C2t = 100*(C2/',num2str(C2,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['C3t = 100*(C3/',num2str(C3,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Bt = 100*(B/',num2str(solution.B(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['wt = 100*(w/',num2str(solution.w(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['rt = 100*(r - ',num2str(solution.R(1)-1,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['tkt = 100*(tauk -',num2str(eco.tauk,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Zt = -ut;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Lt = 100*(Ltot/',num2str(solution.L,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Yt = 100*(Y/',num2str(solution.Y(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['It = 100*(I/',num2str(delta*K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['BoYt = 100*(B/Y*',num2str(solution.Y/((1+eco.tauc)*solution.A-solution.K),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['end;\n\n'] ; fprintf(fid, str);

str = ['%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n'] ; fprintf(fid, str);
%%%%%%%%%%%%%%%% end of the model part


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     Steady-state model       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



str = ['steady_state_model;\n\n'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['a',num2str(h),'        = ',num2str(trunc.a_end_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['at',num2str(h),'       = ',num2str(trunc.a_beg_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    %str = ['c',num2str(h),'        = ',num2str(trunc.c_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['x',num2str(h),'        = ',num2str(x(h),formatSpec),';\n'] ; fprintf(fid, str);
    
end;


str = ['l     = ', num2str(Ramsey.l,formatSpec),';\n'] ; fprintf(fid, str);
str = ['r     = ', num2str(solution.R-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['w     = ', num2str(solution.w,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FK    = ', num2str(1/eco.beta-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FL    = ', num2str(FL,formatSpec),';\n'] ; fprintf(fid, str);
str = ['K     = ', num2str(solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ltot  = ', num2str(solution.L(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['TFP   = ', num2str(1),';\n'] ; fprintf(fid, str);
str = ['A     = ', num2str(solution.A(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['B     = ', num2str(solution.A(1)-solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['u     = ', num2str(0),';\n'] ; fprintf(fid, str);
str = ['tauk  = ', num2str(eco.tauk,formatSpec),';\n'] ; fprintf(fid, str);
str = ['taut   = ', num2str(taut,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ctot  = ', num2str(solution.C,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Y     = ', num2str(solution.Y,formatSpec),';\n'] ; fprintf(fid, str);
str = ['kappa = ', num2str(eco.kappa,formatSpec),';\n'] ; fprintf(fid, str);
str = ['I  = delta*K;\n'] ; fprintf(fid, str);

str = ['ut    = 0;\n'] ; fprintf(fid, str);
str = ['Kt    = 0;\n'] ; fprintf(fid, str);
str = ['Ct    = 0;\n'] ; fprintf(fid, str);
str = ['C1t    = 0;\n'] ; fprintf(fid, str);
str = ['C2t    = 0;\n'] ; fprintf(fid, str);
str = ['C3t    = 0;\n'] ; fprintf(fid, str);
str = ['Bt    = 0;\n'] ; fprintf(fid, str);
str = ['wt    = 0;\n'] ; fprintf(fid, str);
str = ['rt    = 0;\n'] ; fprintf(fid, str);
str = ['Tt    = 0;\n'] ; fprintf(fid, str);
str = ['tkt   = 0;\n'] ; fprintf(fid, str);
str = ['Zt    = 0;\n'] ; fprintf(fid, str);
str = ['Lt    = 0;\n'] ; fprintf(fid, str);
str = ['Yt    = 0;\n'] ; fprintf(fid, str);
str = ['It    = 0;\n'] ; fprintf(fid, str);
str = ['BoYt  = 0;\n'] ; fprintf(fid, str);
str = ['T  = 0;\n'] ; fprintf(fid, str);


str = ['G     = Gss;\n'] ; fprintf(fid, str);
%~ str = ['TLinc = ( G  + (1+FK)*B  - tauk*FK*A - B)/Y ;\n'] ; fprintf(fid, str);
%~ str = ['TKinc = (tauk*FK*A)/Y ;\n'] ; fprintf(fid, str);
%~ str = ['Epsi  = ', num2str(sum(lagMul.psib),formatSpec),';\n'] ; fprintf(fid, str);
%~ str = ['corre = ', num2str(sum(lagMul.psib .* trunc.y0_h)-sum(lagMul.psib),formatSpec),';\n'] ; fprintf(fid, str);

str = ['C1  = ', num2str(C1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['C2  = ', num2str(C2,formatSpec),';\n'] ; fprintf(fid, str);
str = ['C3  = ', num2str(C3,formatSpec),';\n'] ; fprintf(fid, str);


str = ['end;\n\n'] ; fprintf(fid, str);
%%%%%%%%%%% end of the steady state part

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%   Steady-state computation    %%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['resid;\n\n'] ; fprintf(fid, str);
str = ['steady(nocheck);\n\n'] ; fprintf(fid, str);
 
%str = ['check;\n\n'] ; fprintf(fid, str);



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       Aggregate shock        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 str = ['shocks;\n var eps =',num2str(G0,formatSpec),';\nend;\n'] ; fprintf(fid, str);



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%         Stoch simul          %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%~ str = ['stoch_simul (order=1,irf=200) ut Bt tau kappa tauk Yt Ct Kt Lt TLinc TKinc G Y BoYt mut Epsi corre;'] ; fprintf(fid, str);
%~ str = ['stoch_simul (order=1,irf=200) ut Bt tau kappa tauk Yt Ct Kt Lt G Y BoYt mut;'] ; fprintf(fid, str);


str = ['stoch_simul (order=1,irf=500,periods=1000) ut Bt taut kappa tauk Yt Ct It Lt C1t C2t C3t T;'] ; fprintf(fid, str);

isOctave && fflush(fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        Launch Dynare         %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    command = ['dynare ',file, ' noclearall'];
    eval(command);

name = ['fig_trunc',num2str(rho_u),'.mat'];
%save (name,'ut_eps','Bt_eps','tau_eps','kappa_eps','tk_eps','Yt_eps','Ct_eps','Kt_eps','Lt_eps','TLinc_eps','TKinc_eps','G_eps','Y_eps','BoYt_eps','mut_eps','oo_')


save (name,'C1','C2','C3','ut_eps','Bt_eps','Yt_eps','It_eps','Lt_eps','Ct_eps','C1t_eps','C2t_eps','C3t_eps','kappa_eps','oo_')




